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ABSTRACT 

in 

In a further development of a deterministic planet-formation model (Ida & 
Lin 2004), we consider the effect of type-I migration of protoplanetary embryos 
due to their tidal interaction with their nascent disks. During the early embedded 
phase of protostellar disks, although embryos rapidly emerge in regions interior 
to the ice line, uninhibited type-I migration leads to their efficient self-clearing. 
But, embryos continue to form from residual planetesimals at increasingly large 
radii, repeatedly migrate inward, and provide a main channel of heavy element 
accretion onto their host stars. During the advanced stages of disk evolution (a 
few Myr), the gas surface density declines to values comparable to or smaller 
than that of the minimum mass nebula (MMSN) model and type-I migration is 
no longer an effective disruption mechanism for mars-mass embryos. Over wide 
ranges of initial disk surface densities and type-I migration efficiency, the surviv- 
ing population of embryos interior to the ice line has a total mass several times 
that of the Earth. With this reservoir, there is an adequate inventory of residual 
embryos to subsequently assemble into rocky planets similar to those around the 
Sun. But, the onset of efficient gas accretion requires the emergence and reten- 
tion of cores, more massive than a few M e , prior to the severe depletion of the 
disk gas. The formation probability of gas giant planets and hence the predicted 
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mass and semimajor axis distributions of extrasolar gas giants are sensitively 
determined by the strength of type-I migration. We suggest that the observed 
fraction of solar-type stars with gas giant planets can be reproduced only if the 
actual type-I migration time scale is an order of magnitude longer than that de- 
duced from linear theories. We also show that the introduction of such slower 
type-I migration rate makes simulated planetary mass-period distribution more 
consistent with observation. 

Subject headings: extrasolar planets - planetary systems: formation - solar sys- 
tem: formation 



1. Introduction 



In the radial velocity survey of nearby stars, more than 200 extrasolar planets with mass 
M p comparable to that of Jupiter Mj have been discovered. Current statistics indicates that 
i]j > 7% of all nearby F, G, and K dwar fs on various sear c h target lists have gas giant 



planets with semimajor axis < 5AU (e.g., iMarcy et al.l 120051 ; iMayor et al.ll2005l ). But the 



extrapolation of existing data suggests a h igher fraction of solar-type stars may have longer- 
period gas giant planets (jCummingi 12004 ) . 



In the previous papers of this series (llda fc Lin ll2004aU bl. 120051 hereafter referred to as 
papers I, II, and III), we carried out several sets of simulations with a numerical prescription, 
which is based on the sequential accretion scenario. We used these results to explain the 
observed mass and orbital distributions of extrasolar planets and to constrain the intrinsic 
physics of planet formation through detailed comparisons with the observed data. 

In these previous investigation, we have taken into account the effect of gas depletion 
by assuming that the gas surface density S s declines everywhere exponentially over a char- 
acteristic timescale Td ep — 1-10 Myr. With an a prescription for the turbulent viscosity, this 
assumed global evol ution of S fl can lead t o declining accretion rates which are consistent 
with those observed (IHartmann et al.lll998l ). T he magnitude of t^ p sets a strong constraint 



on the gas giant planet building efficiency (e.g.. llda fc Lin Il2004al ). 



We have also assumed that the dust surface density (S^) in these disks is preserved from 
the initial value we have adopted. The dust grain s' mass inferred from the m m observation of 
continuum radiation from protostellar disks (e.g.. lBeckwith fc Sargentlll996l ) appears to have 
a wide dispersion cen tered around the value of the minimum-mass solar nebula (MMSN) 
model (IHayashil Il98ll ). We adopted a similar distribution for centered around that of 
MMSN. This assumption, though it greatly simplified our treatment, is less well justified. In 
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the scaling of with that of MMSN, we have inherited the assumption that all the heavy 
elements were locally retained during the epoch of planet formation. 

There are two potential avenues for the depletion of heavy elements: hydrodynamic gas 
drag of small dust grains and type-I migration of cores due to the tidal interaction with their 
nascent disks. Here, " cores" mean the protoplanetary e mbryos which formed as a result of 
runaway/oligarchic growth of accretion of planetesimals ( IKokubo fc Idalll998l . 120021 ) . These 
cores are generally not sufficiently massive to initiate runaway gas accretion. Without much 
additional growth before the disk gas is severely depleted, these cores would become either 
rocky planets or ice giants. 

In the present paper, we consider the dominant processes after the formation of planetary 
embryos which are hydrodynamically decoupled from disk gas motions. Prior to the phase of 
rapid gas accretion, these embedded cores a nd their surrounding gas exchange angular mo- 
ment um as they engage in tidal interaction (IGoldreich fc Tremaind Il980l ; iLin fc Papaloizou 
19791 ). Analytic studies suggest that isolated cores lose angular momentum to the disk ex- 
terior to their orbits faster than that they gain from the disk interior to their orbits. This 
torque imbalance leads to a "type-I" migration. From the linear analysis, the characteristic 
orbital deca y time scale of Earth-mass core s at several All's in a MMNM is much shorter 
than 1 Myr (1 Wardlll986l ; iTanaka et al.ll2002l ). Accordingly, embryos formed well outside the 
ice line may migrate to the present locations of the Earth and Venus. But, the mostly re- 
fractory compositions of the terrestrial planets in the present-day solar system suggest that 
they probably did not migrate to their present location from regions well beyond the ice line. 

Today, there remains considerable uncertainties in the efficiency of this process (see 
§2.2). Nevertheless, their formation could have been preceded by a much larger population 
of cores which did migrate into the Sun. In view of this uncertainty, we carry out a parametric 
study on the decline of the cores' accretion rate associated with a T, d reduction due to the 
type-I migration of cores (see below). This effect has been neglected in our previous papers. 

Heavy elements are also depleted from the disks through "type-II" migration. Due to 
runaway gas accretion onto a core, the growing gas giant planet eventually acquires sufficient 
mass to open a gap in the disk and lock itself into a "ty pe-II" migration with the viscous 
evolution of th e disk gas (Lin fc Papaloizoul Il985l . 119931 ). In the context of solar system 
formation. ILin I (119861 . 119951 ) has speculated that, when the solar nebula was massive, several 
gas giants may have formed, undergone type-II migration, and eventually merged with the 
Sun. When the disk's S g decreases below that of the MMSN model, the migration has 
stalled. According to this scenario, Jupiter and Saturn could be the last survivors of several 
generations of protoplanets (also see iPapaloizou and Terqueml 120061 ). although it may not be 
easy to form terrestrial planets after preceding gas giants have passed through the inner disk 
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region (lArmitage 1 120031 ) . A resurgence of this migration and survival scenario (jLin et al. 
19961 ) have been stimu l ated by the obser vational discoveries of close-in gas giant planets 



'e.g.. lMarcy et al.ll2005l ; iMayor et al.ll2005l ). But, a majority of the known extrasolar planets 
have periods much longer than a few days. We showed in Papers I and II that their observed 
logarithmic period distribution can be used to infer comparable timescales for giant planets' 
migration and disk depletion. However, based on the relative abundance between planets 
with a < 0.06 AU and those between 0.2-2AU's from their host stars, we have to assume 
that up to 90% of the gas giants which migrated to the proximity of their host stars may 
have perished. 

Here, we focus on the repeated clearing due to type-I migration. We suggest that 
this self-regulated clearing mechanism limits and determines the amount of heavy elements 
retained by the cores and residual planetesimals (also see discussion on metallicity homo- 
geneity of open clusters in §4.1). It also results in late formation of gas giants and the 
marginal probability of gas giants' formation through a series of failed attempts. 



The scenario we consider here is similar to the hypothesis proposed by ICanup fc Ward 



(120061 ) in the context of satellite formation around Jovian planets. They suggested that 
the total mass of the surviving satellites is self-regulated by their type-I migration through 
their nascent cir c um-planet ary disk which is continuously replenished. In the context of 
planet formation, iMcNeil et al.l (120051 ) carried out N-body simulations of terrestrial planets' 
accretion and growth including the effect of type-I migration. They found that it is possible 
to retain sufficient amount of solid materials (planetesimals and embryos) to assemble Earth- 
mass planets within a disk lifetime ~ 10 6 years provided the type-I migration rate is slightly 
slower than that predicted by the linear theory. In these simulations, the initial disk mass 
in terrestrial planet regions is assumed to b e 3-4 times o f that of MMSN. Based on an 
analytic approximation for type-I migration, iDaisaka et al.l (120061 ) simulated the evolution 
of T, d under various disk conditions. Their results indicate that the depletion region expands 
outwards, starting from the disk's inner boundary. They also found that in disks with same 
dust-to-gas ratio, the asymptotic surface density of the retained embryos decreases with the 
initial value of (and £ 9 ) because type-I migration is faster in more massive disks. 

Here we consider a much large range of initial conditions. With a comprehensive nu- 
merical prescription, we consider the possibility of multiple generation of embryo formation. 
We show that in relatively massive disks although type-I migration is more effective for indi- 
vidual planetesimals and embryos, it does not necessarily lead to more rapid depletion of the 
residual population. In these disks, migration starts with smaller individual masses but more 
generations of embryo form and parish before the initial supply of planet-building blocks is 
severely depleted. Our results indicate that the total retainable mass of heavy elements only 
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weakly depends on the initial disk mass. 



Although type-I migration places a mass limit on the retainable embryos, it does not 
quench the formation of Earth-like terrestrial planet s because they can b e assembled from re- 
tainable low-mass embryos after the gas is depleted (IDaisaka et al.ll2006l ). Gas giant planets, 
however, must form in gas-rich disks which appears to be depleted on the time scale Td ep ~ 
a few Myr. In addition, they must acquire ~ 10M e cores prior to the onset of efficient gas 
accretion. Although they can form rapidly outside the ice lines on ma ssive disks, early gen- 
eratio n of such massive cores quickly migrate into their host stars. But, iThommes fc Murray 
( 120061 ) showed that cores formed after the disk gas has been severely depleted may with- 
stand the disruption by the declining type-I migration. In §3 of thi s pape r, we present results 
which are consistent with that obtained by IThommes fc Murray! (120061 ). We also incorpo- 
rate gas accretion onto cores and the effect of type-II migration and find that under some 
circumstances, there is adequate residual gas to promote the efficient gas accretion and the 
formation of gas giant planets. This "late formation" scen ario is conceptually co nsistent 
with that inferred from the noble gas enrichment in Jupiter (IGuillot fc Huesoll2006l ). Based 
on this model, we are able to reproduce the observed mass-period distribution of the known 
gas giant planets around solar type stars provided the efficiency of type-I migration is at 
least an order of magnitude slower than that derived from the linear theory (see §3). Our 
conclusion is consistent with the studies by ? who also found that a necessary condition for 
the formation of Jupiter and Saturn is a 30-times reduction in the type-I migration rate. 

In metal-poor disks, the formation of critical mass cores requires a relatively large gas 
surface density. In these disks, type-I migration is more effective in clearing the cores prior 
to the onset of gas accretion. Consequently, the formation of gas giant planets is suppressed. 
Assuming the disks' initial metallicity is identical to that of their host stars, this effect 
sharpens the dependence of the predicted formation efficienc y of gas giants on the st ellar 
metallicity which is well established in the observational data (IFischer fc Valenti II2005I ). 



Since type-I and II migrations are the essential processes which determine the properties 
of emerging planets, we briefly recapitulate, in §2, their basic physical principles. We also in- 
corporate a quantitative prescription of these processes in our existing comprehensive model 
for planet formation. The model for the new addition, type-I migration and the decrease in 
planetesimal surface density due to the migration, is explained in detail. In §3, we carry out 
a systematic study on how cores' migration may affect the formation of terrestrial planets 
and gas giant planets. We show that type-I migration delays formation of gas giants, and 
that it leads to the mass and semimajor axis (M p — a) distribution consistent with that of 
the known gas giants around G dwarfs. We also study dependence of the the distribution on 
planet formation parameters, metallicity. The introduction of relatively slow type-I migra- 
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tion results in the metallicity dependence that may be consistent with observation. Finally, 
we summarize our results and discuss their implications in §4. 

2. Planet formation and migration model 

In this paper, we simulate the M p — a distribution of extrasolar planets, taking into 
account the effects of cores' type-I migration as well as type-II migration. In Papers I-III, 
we outlined in detail a quantitative prescription which we used to model the evolution of 
planetesimals, gas accretion onto protoplanetary cores, the termination of gas giant planet 
growth, and type-II migration. In this section we briefly recapitulate the dependence of both 
types of migration on the background disk properties. 

In this paper, we newly incorporate the effects of type-I migration and examine the 
impact of type-I migration on the formation of terrestrial planets and the asymptotic prop- 
erties of emerging gas giant planets. In light of the theoretical uncertainties, we introduce 
a prescription which captures the main determining factors of type-I migration and enables 
us to consider a wide range in the magnitude of its efficiency factor. We here consider only 
the case of M* = 1M to focus on the effects of type-I migration around solar-type stars, 
although we retain the dependences on stellar mass M* and its luminosity L*. 

Note that inclusion of type-I migration and the minor changes in truncation conditions 
of growth of gas giants and the disk model that are described below do not change the 
conclusions that have been derived in previous papers: (i) the existence of "planet desert" (a 
lack of intermediate mass planets at < a few AU) in Paper I, (ii) the metallicity dependence 
(a fraction of stars with giant planets increases with the metallicity) in Paper II, and (iii) 
the stellar mass dependence (giant planets are much less abundant around lower-mass stars) 
in Paper III. The results with inclusion of type-I migration corresponding to (i) and (ii) are 
shown in Figures [5] and [6], respectively. The stellar mass dependence with type-I migration 
is briefly commented on in §3.2.3 and will be discussed in detail in a separate paper. 

2.1. Parameterized disk models 

The main objectives of this series of papers is to examine the statistical properties 
of emerging planets under a variety of disk environments rather than to study individual 
processes which regulate the disk structure. For computational convenience, we introduced, 
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in Papers I- 1 1 1 , multiplicative factors (fd and f g ) to scale S 5 and such that 



Z d ^^(r/lOAU)"*, 

S g = £ ff , 10 / fl (r/10AU)-*, 1 J 

where normalization factors £^10 = 0.32g/cm 2 and S 9l io = 75g/cm 2 correspond to 1.4 times 
of S g and at 10AU of the MMSN model, and the step function ?7i ce = 1 insi de the ice line 
at ajp e (eq. ^\ 
1994T )]. 



and 4.2 for r > a ice [the latter can be slightly smaller (~ 3.0) (jPollack et al. 



We show below that the disk metallicity [Fe/H]^ is an important parameter which regu- 
lates the survival of protoplanetary cores during their type-I migration. Dependence of disk 
metallicity is attributed to distribution of fd : o = / g> olO^ Fc / H ^. Solar metallicity corresponds 
to [Fe/H]d = and fd,o = f g ,o- The main advantage of these parameterized disk structure 
models is that they enable us to efficiently simulate and to examine the dominant dependence 
of planet formation and dynamical evolution on the disk structure. 

In self-consistent treatment of accretion flow, the disk temp erature is determined by a n 
equilibrium between the viscous dissipati on and heat transport (IShakura fc Sunvaevi 1973). 
We neglect the detailed energy balance ( Chiang fc Goldreich 1997 ; Garaud fc Lin 12007 ) 
an d adopt the equ ilibrium temperature distribution in highly optically thin disks prescribed 
by lHayashil (1198ll ) such that 

1/4 

T = 280 ( - — r~rr ) I ^ 1 K. (2) 



1AU 



In this simple prescription, we set the ice line to be that deter mined by an equilibrium 
temperature (eq. [2]) in optically thin disk regions ( Hayashil 1981 ). 



2.7(L,/L ) 1 / 2 AU, 



(3) 



where L* and L Q are the stellar and solar lumino sity. The magnitude of a\ re , may be modified 
by th e local viscous dissipat ion (ILecar et al.ll2006l ) and stellar irradiation (j Chiang fc Goldreich 
19971 ; iGaraud fc Lin 1 120071 ). These effects do not greatly modify the disk structure during 
the late evolutionary stages and they will be incorporated in subsequent papers. 



2.1.1. Evolution of f g 



1985|), 



The surface density of the gas can be determined by a diffusion equation (ILin fc Papaloizou 



<9£ 9 _ 1 d 
dt r dr 



3r^|-(EX /a ) 



0. 



(4) 
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where we neglected the sink terms due to photoevaporation and accretion onto the cores. 
For computational conveni ence, we adopt the standard constant-a prescription for viscosity 
( IShakura &: Sunyaevlll973l ). 

v = ac s H, (5) 

where c s , H = v^2c s /f2K, and are the sound speed, disk scale height, and the Keplerian 
angular frequency. Although the value of a is not clear, a ~ 1 0~ 3 may be consistent with 
the observational data of accretion r ates of T Tauri d isks (e.g., lHartmann et al.lll998l ) and 
the results of MHD simulations (e.g.. lSano et al.ll2004l ). We here use a = 10~ 3 as a nominal 
value. 

With the a prescription and eq. (121), v oc r and eq. (j4j) reduced to a linear partial differen- 
tial eq uation for which self-similarity solutions have been presented by lLynden-Bell fc Pringle 
(11974] ) . The numerical solution for disk-evolution equation [i.e., eq. (@J with a = 10~ 3 ], start- 
ing with S 9 oc r~ L5 for r < r m and an exponential cutoff at r m = 10AU, is illustrated in 
Figure [Th- 

After a brief initial transition, the numerical solution quickly approaches to £ 9 oc 
r -1 with an asymptoti c expo nential cut-off, which is the self-similar solution obtained by 
Lvnden-Bell &: Pringld (119741 ) for a linear viscosity prescription and by iLin &: Bodenheimer 



for an a model. As shown in Figuredb, S 9 in the region r < r m decreases as uniformly 
independent of r. In the self-similar solution, S 9 at r < r m decays as S g oc (t/r dcp + 1)~ 3 / 2 
where 



cp 



3 x 10 fc 



a 



Mrm) V10- 3 
In the above, we used eqs. ([2]) and 



1 ( Tm \ ( M * 
VlOOAUy \M Q 



-1/2 



yrs. 



(6) 



This self-similar nature of the solution is preserved if we assume a spatially uniform 
exponential depletion of the disk (i.e. S 9 oc exp(— t/rd ep ) with a r-independent Td ep ), as in 
Papers I-III, although in the actual self-similar solution, S 9 oc t -3 / 2 for t > t^ cp . If the effect 
of photoevaporation is taken into account, T, g decays rapidly after it is significantly depleted, 
so that the exponential decay could be more appropriate. Thus, we adopt 



f g = / g ,oexp(-t/r d , 



epji 



(7) 



in the disk model in eq. ([T]). Even if we fix the a value, Td ep has uncertainty, because we 
do not have enough knowledge about r m , and r m should have some dispersion. We use Td ep 
as a parameter and set it to be in a range of 10 6 -10 7 yrs, which may be consistent with 
observation. Corresponding to the self-similar solution, we consider here disk models with 
Tjg oc r _1 (q g = 1) in addition to the case with E 9 oc r~ 3//2 (q g = 3/2) that was assumed in 
Papers I-III corresponding to the MMSN model. 
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2.1.2. Evolution of f d 



We start each model with E^ oc r -3 / 2 (qd = 3/2). In principle, the evolution of E^ 
should be treated independently from E 5 even though small dust particles are thermally and 
dynamically coupled to the gas. Physically, the magnitude of E^ is determined by the grains' 
growth rate and the planetesimals' retention efficiency In view of the large uncertainties 
in these processes, we adopt the radial gradient of conventional MMSN model, qd = 3/2. 
Different radial gradient between E^ and E 9 could be produced by inward migration of dust 
grains due t o gas drag, which tends to make inner disk m ore metal-rich while outer disk 
metal-poor (IStepinski fc Valageaslll997l ; iKornet et al.ll200ll ). 



In Papers I-III, fd was assumed to be constant with time until core mass reaches iso- 
lation mass. In the present paper, we take into account planetesimal clearing due to the 
cores' accretion. Since we also include type-I migration of the cores, fd at a given location 
(semimajor axis a) continuously decreases with time as planetesimals are accreted by cores 
which in term undergo orbital decay. Note that in the present paper, a is identified as r, 
since we neglect evolution of orbital eccentricities. 



The full width of feeding zone of a core is given by (IKokubo fc Idalll998l . 120021 ) 



Aa r 



10n 
10 



2M C 
AT. 



1/3 



where rn is two-body Hill radius. Following Papers I-III, we take into account of the expan- 
sion of fe eding zones due to coll i sions among the isolated cores after significant depletion of 
disk gas (IKominami fc Idall2002l ; IZhou et al.ll2007l ). although this effect influences the results 
only slightly (it is important in inner regions in which cores are significantly depleted by 
type-I migration). 

An increase in the cores' mass is uniformly subtracted from the mass in its feeding 
zone, keeping qd locally. When the cores migrate out from the feeding zone, the reduction 
of fd at the initial location of the cores is stopped. But, the cores can continue to accrete 
planetesimals along their migration path. In any given system, when a core reaches 0.03AU 
(occasionally, it can attain sufficient mass to accrete gas and become a gas giant during the 
course of its migration), a next-generation core is launched at the pre-migration radius with 
100 smaller mass. (Since core growth is faster during earlier stage of disk evolution (eq. [9]), 
the choice of initial mass of the next-generation core does not affect the total core accretion 
timescale.) Thereafter the growth of the core and the depletion of the planetesimals in its 
feeding zone resume. 



This process repeats until the residual planetesimals are depleted. Along the cores' 
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migrating paths, they accrete residual planetesimals (until they reach a < ad ep ,mig given by 
eq. [27]; see below) but cannot reduce fd as fast as the cores formed in situ. Therefore we 
neglect the evolution of fd due to the planetesimal accretion by the migrating cores. The 
clearing of the residual planetesimal disk by migrating gas giants, is also neglected because 
it rarely occurs and is unlikely to affect the overall distribution of the emerging terrestrial 
planets (see §3.2.2). The overall evolution of due to core accretion and migration is 
described in 53.1. 



2.2. Core growth and type-I migration 

We briefly recapitulate a comprehensive analysis on the growth of planetesimals and 
cores. Readers can find a systematic derivation of these results and the appropriate references 
from Paper I. During their growth through cohesive collisions, cores' mass accretion rate (of 
planetesimals) at any location a and time t is described by (Paper I) 

dMjdt = M c /r c>acc ; 

r c , acc = 2.2 x lO^/,- 1 /,- 275 ^^^-^) ( _ | _ ) -/io + (, d -3/ 2 ) + ( 2 / 5 ) fe -3/ 2 ) ^)" 1/6 yrS) 

(9) 

where fd and f g change with time. In the derivation of the above expression, we have 
adopted the mass of the typical field planetesimals to be m = 10 20 g. The numerical factor 
1q(2/5)(3/2— c/ 3 ) comes f rom the fact that we scale surface densities at 10AU in eq. ([1]). 

Based on the results of previous numerical simulations (see references in Paper I), 
we assume that runa way growth phase is quickly transferred to oligarchic growth phase 



( iKokubo fc Idalll998l ). During the oligarchic growth phase, the cores' accretion timescale is 

1 /3 

dominated by their late-stage growth (r Cjacc oc M c ) so that dependence on their initial mass 
is negligible. In eq. Q), we adopt M c (0) = m = 10 20 g. When f g < 10~ 3 , we use gas-free 
accretion rate rather than that in Eq. (Q (Paper I), 

7 i , i ( a \ 3+fe-3/2) / MA ~ 1/2 , . 

Vc = 2xl0 7 ^ __ -L) yrs . (10) 



1AU/ \M, 



1 G 



In Papers I-III, we artificially terminate the cores' accretion at M c = M Ci i SO , where the 
isolation mass is given by 

nlR 3/2,3/2/ a x3/4-(3/2)( gd -3/2) / M ,y 1/2 

In the present paper, we do not need to adopt an artificial termination to the core growth 
because when M c approaches isolation M ciso , in situ core accretion is automatically slowed 
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down by the reduction of fd- However, in outer regions, ejection of planetesimals by large 
protoplanets limits their accretion at M csca (scattering limit; Paper I). We take this effect 
into account by adopting an artificial termination for the cores' accretion of planetesimals 
when they attain a mass M CiSca ~ 1.4 x lO 3 (a/lAU)~ 3 / 2 (M*/M )M e . 

Imbalance in the tidal torques fr om outer and i nner d isks causes type-I migration of a 
core. Through 3D linear calculation, iTanaka et al.l (120021 ) derived the time scale of type-I 
migration, 



migl 



a d 2.728+1.082g 9 

5 x 10 4 x 1Q( 3 / 2 -^) 



Olfg 



m» M t n-i 

Mr, o^Sl K 



Ms 



I a \1a 

Uau/ 



M*.\ 

M Q J 



3/2 



(12) 



yrs. 



where we used eq. (T5]) and q g = 1-1.5. 



In eq. (T121 . we introduce an scaling factor C\ to allow for the retar dation of type-I 



migration due to non-linear effects. In the expression of ITanaka et c 



simulations have been carried out (see e.g.. |Papaloizou and Terquem 



1|2002|), d = 1. Many 



20061 ) with different nu- 



merical methods and resolutions, protoplanetary potentials and orbits, and disk structures. 
They produced a considerable range of the timescale of type-I migration (r migl ). Retarda- 
tion processes f or this type-I migrat ion include variation in the surfac e density and temper- 



ature gradient (|Masset et al.l l2006bl ) . intrinsic turbulence in the disk (|Laughrin et al 



Nelson fc Papaloizoul 120041 ). self-induced unstable flow flKoller fc Li 



200 



2004 



Li et al.1 120051 ) 



and non-linear radiative and hydrodynamic feedbacks (IMasset et al.l l2006af ). Under some 
circumstances, C\ is reduced to be < 0.1 (Dobbs-Dixon et al. 2007, in preparation). In light 
of these uncertainties, we adopt C\ as a parameter and mainly examine the slower migration 
cases with C\<\. 

Type-I migration can occur before core mass reaches M c iso . In this case, (equiv- 
alently, fd) decreases but does not completely vanish. From the decreased S^, new cores 
accrete with reduced growth rates and isolation masses. The cycle of core growth, type-I 
migration, and reduction of con tinues until decr eases to the levels from which only 
cores small enough not to migrate (jDaisaka et al. 2006 ). We present a detailed discussion 
on the self-regulation mechanism in later sections. 

Under some condit ions, the migrating cores can capture an d accumulate planetesi- 
mals along th e ir pat hs (IWard fc Hahnl Il995l ; ITanaka fc Idalll999l ). N-body simulation by 
Daisaka et al.l (120061 ). however, showed that the trapping of planetesimals by the cores is 
tentative and it does not significantly reduce their accretion rates. In our simulations, we 
use the accretion rate for migrating cores that is the same as the rate for non-migrating cores 



(eq. [9]). For cores in the systems with q g — 1 and qd = 1.5, 

dMjdt = M c /r c , acc ; 

r c , acc = 3.5 x 10^-7," V /5 (mj) 5/2 (ife)~ V Vs, 
da/dt = a/r mig i ; 

r migl * 1.6 x lO^fV- 1 (1^) (ft)' 72 yrs 

The magnitude of fd is consistently decreased by the increase in M c (see §3.ip . Since 
the accretion rate is determined by the instantaneous local value of fd, it limits the mass of 
the migrating cores in a gaseous medium. Prior to the severe depletion of the disk gas, we 
quench the cores' accretion of planetesimals at a < ad ep ,mig (eq. [27]) on the basis that there 
is an inadequate supply of residual planetesimals in these locations to significantly add to 
their masses. But, during the gas depletion, we assume f g decays exponentially. After disk 
gas is significantly depleted (f g < 1CT 3 ), Eq. (TTDT) is used for r cacc . We set also da/dt is set 
to be zero at disk inner edge (~ 0.03AU). 



(13) 
(14) 



2.3. Formation of gas giant planets 



Prescriptions for formation of gas giant planets are the same as those used in Paper 
I-III, although the cores' rate of planetesimal accretion is revised from Paper I-III by the 
incorporation of type-I migration and planetesimal depletion. After the formation of these 
gas giants, we assume all residual planetesimals in the gap is cleared as consequence of 
dynamical instabilities. The se planetesimals are either accreted by the gas giants or scattered 
elsewhere ( IZhou et al.l 120071 ). We neglect the emergence of second-generation cores close to 
the orbit of gas giants. 

In principle, cores with mass much less than that of the Earth can accrete gas. But, 
unless the heat released during gas and planetesimal accretion is diffused and radiated away, 
quasi thermal and hydrodynamic equilibrium would be established to prevent further flow 
onto the cores. Around low-mass cores, the temperature and density of the envelope are low 
so that heat cannot be easily redistributed through their envelopes. But as the cores grow 
through planetesimal bombardment beyond a mass 

0.25 

M c , hydro ~ 10 ( TT^nrr^z I M e, (15) 



lO^Me/yr^ 

both the radiative and c onnective transpor t of heat become efficient to allow their envelope 
to contract dynamically (jlkoma et al.ll2000l ). In the above equation, we neglected the depen- 
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dence on the opacity in the envelope (see Paper I). In regions where the c ores have already 
attained isolation, their planetesimal- accretion rate M c is much diminished ( Zhou et all2007 ) 



and M c hydro can be comparable to an Earth mass M e . But, gas accretion also releases energy 
and its rate is still regulated by the efficiency of radiative transfer in the envelope such that 



dM p 



TKH 



(16) 



where M p is the planet mass including gas envelope. In Paper I, we approximated the 
Kelvin-Helmholtz contraction timescale tkh of the envelope with 



t kh ?t 10 



A.-1 



-k2 



yrs. 



;i7) 



In order to take into account the uncertainties associated with planetesimal bombardment, 
dust sedimentation and opacity in the envelope, we adopt a range of values kl = 8-10 and 
k2 = 3-4 (see Papers I and II). Here we use k2 = 3 and treat kl as a parameter. 

Gas accretion onto the core is quenched when the disk is depleted either locally or 
globally. We assume that gas accretion is terminated if either the thermal condition or 
global depletion condition is satisfied. A (partial) gap is formed when the rate of tidally 
induced angular momentu m exchange by the plane t with the disk exceeds that of the disk's 



intrinsic viscous transport (ILin &: Papaloizou 



1983) (the viscous condition) 



M p > M g 



30 



a 



10" 3 / V1AU 



1/2 



Planets with M p > M gjVis have sufficient mass to induce the partial clearing of the disk 
near their orbit. The tidal torque on either side of the gap becomes sufficiently strong to 
induce the planets to adjust their positions within the gap. This feedback process leads to 
a transition from type-I to type-II migration. Therefore, we adopt the viscous condition for 
the onset of type-II migration. In Papers I- 1 1 1 , we adopted a = 10~ 4 as a nominal value. 
In order to match the simulated M p — a distribution to the observed data, eq. (ITS]) was 
arbitrarily multiplied by a scaling factor A v = 10. 

The previously adopte d value of a is smalle r than that inferred from the models of 
protostellar disk evolution (IHartmann et al.lll998l ). In this paper, we use a = 10 -3 without 
imposing the scaling factor. As a result, the condition for the onset of type-II migration is 
the same as that in Papers I-III. 

A clear gap would be formed and gas accretion would be terminated wh en the planet's 
Hill radius becomes larger than disk scale height (ILin fc Papaloizou! 119851 ) (the thermal 
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3/4 / L* \ 3/8 / M* N _1/2 



condition) and it is given by (Paper I) 

M r > M g , th = 0.95 x 10= j j (Jg j M e . (19) 

While the thermal condition rn > 1.5if was used in Papers I- 1 1 1 , we here used the condition 
?"h > 2H, in order to make clear that the asymptotic mass of gas giants is determined by 
global depletion of disk gas rather than the local thermal condition in the cases with inclusion 
of type-I migration in which gas giants are formed in relatively late stages. (It can also reflect 
the effect of gas flow into the gap mentioned below.) 

Some numerical simulations indicate that a residual amount of gas may continue to 



flow i nto the gap after both the viscous and thermal c onditions are satisfied (lD'Angelo et la. 



20031 : Kiev fc Dirksen 



2006 



Tanigawa &: Ikomal 120071 ). However, recent numerical simula- 



tions (IDobbs-Dixon et al.ll2007l ) also show that the azimuthal accretion flow from the coro- 
tation regions onto the planet is effectively quenched despite a diminishing flux of gas into 
the gap. Since the thermal condition usually requires larger M p than the viscous condition 
(M gjt h > M giV ; s ), we adopt it as the criterion for the termination of gas accretion in the 
determination of the asymptotic mass of gas giant planets. But, we take into a residual 
amount of gas which may leak through the gap and provide an effective avenue for angu- 
lar transfer between the inner and outer regions of the disk via the gas giant's corotation 
resonance (§2.4). 

Gas accretion may also be limited by the diminishing amount of residual gas in the 
entire disk even for planets with M p < M P)t ^. In Papers I-III, we assumed that the maximum 
available mass is determined simply by M gnoiso ~ 7ra 2 S 9 . However, because the global limit 
plays a more important role when type-I migration is incorporated, we use a more appropriate 
condition M g . 

noiso ~ Jq 2,TTQjYjgdct, (We neglect disk gas inflow from far outer regions). For 

% = h 

M g>noiso ~ 3.5 x 10 2 / g ,oexp (j^u)^ M ©- ( 20 ) 

When M ginoiso diminishes below M p or when M p exceeds M gjt h, gas accretion is terminated. 



2.4. Type-II migration 

In our previous analysis in Papers I-III, we adopted a simply analytic prescription for 
type-II migration: (i) before disk gas mass decays to the value comparable to the planet 
mass, the planet migrates with (unperturbed) disk accretion and (ii) when disk gas mass is 
comparable to the planet mass, a fraction (C2 ~ 0.1) of the total (viscous plus advective) 
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angular momentum flux transported by the disk gas (which is assumed to be independent 
of the disk radius) is utilized by the planet in its orbital evolution. 

If the planet's tidal torque can severely clear a gap in the vicinity of its orbit, £ 9 in inner 
disk region would decrease faster than that in the outer region. The full torque asymmetry 
leads to G% ~ 1 in case (ii). But, when the truncation condition is marginally satisfied, 
the disk interior to the planet's orbit have a total mass > M p such that it can effectively 
replenish the angular momentum lost by p lanet to the outer di sk region. There may also 
be a leakage of gas through the gap region (jD'Angelo et la.ll2003l ) which would suppress the 
degree of torque asymmetry. The protoplanet's corotation resonance may drive an effective 
angular momentum transfer across the two disk regions separated by the protoplanet. 

Protoplanets are formed in disk regions where the midplane is inactive to magneto- 
rotational instabilities. It is possible that the modest accretion rate onto the host stars flow 
through this region yia an active la yer which is exposed to external ionizing photons and 
cosmic ray particles (IGammie Ill996l ). Different values of effective ct's may contribute to the 
mass flow through the disk and the planet-disk interaction, especially if there is some leakage 
across the gap region. All of these possible scenarios can lead to C2 <C 1. 

The uninterrupted replenishment of gas into the corotation region may also mai ntain a 
finite vortensity gradient and an unsaturated corotation torque fjMasset et al.ll2006af ) which 
may reduce the efficiency of type-II migration from disk gas accretion even in case (i) 
( jCrida &: MorbideUil 120071) . However, the results of another set of 2D numerical hydro- 
dynamic simulations (jD'Angelo et al.l 120061 ) essentially reproduces the the ID simulation 
( )Lin fc Papaloizoul Il986l ) in which the contribution from the corotation resonance is ne- 
glected. 

In order to take into account of these uncertainties, we reduce here the migration rate 
by a factor C 2 for cases (ii) as in Papers I-III, 



^~mig2,ii 



(l/2)M p Q K (a)a 2 



C 2 37rS 9 (r m )r^z/ m f2 K ( 



5xio 5 rM^ 
10 



-1 



Mi 



a \ 1/2 / M* 



M, 



1AU 







-1/2 



(21) 



yrs. 



When T mig2ji i is shorter than the migration timescale for (i), we use the latter timescale, 



7~mig2,i 



a 1 (3/2)i/ /a 



0.7 x 10 5 



a 



-1/2 



(22) 



yrs. 



.10-V VlATjy \M Q , 
Due to more accurate estimations, the numerical factors in r mig2i i and r mi g 2jii slightly differ 
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from those in Papers I-III for the same values of C 2 and a. However, there are still uncertain- 
ties in the formula for T mig2i ii. Furthermore, it is not easy to theoretically evaluate the value 
of C 2 as well as a. Hence, we vary the magnitude of C 2 and compare the simulated results 
with the observed data to obtain a calibration. For a ~ 10~ 3 , in the limit that C 2 ~ 1, most 
of gas giants are removed from the regions beyond 1AU, which is inconsistent with observed 
data of extrasolar planets. As shown in Paper II, in order to reproduce M p -a distribution 
of observed extrasolar planets, disk depletion timescale r dep must be ~ T mig2 at a few AU. In 
the present paper, we mostly adopt C 2 a = 10~ 4 . For a ~ 10~ 3 , it corresponds to C 2 ~ 1/10. 
As a result, r mig2) ii we adopt here is almost identical to that used in Papers I-III. 

3. Effects of type-I migration 

In our previous simulations in Papers I-III, the effect of type-II migration was included 
but that of type-I migration was neglected. In this section, we investigate the effects of 
type-I migration on the planets' M p -a distribution. 

3.1. Surface density evolution due to type-I migration 

In order to analyze the asymptotic M p —a distribution from the Monte Carlo simulations, 
we first present the results on the reduction due to type-I migration. The relevant time 
scale here is a function of C%, f g and M p (eq. [12] )• The magnitude of the cores' M p is a 
function of fa, f g and t. 

Since type-I and II migrations involve the tidal interaction of embedded cores with the 
disk gas, the dynamical clearing of the residual planetesimals is suppressed after the gas 
is severely depleted (except for outer regions in which ejection by massive embryos can be 
efficient). Since t ~ T^ ep is a critical stage for the build-up and retention of the cores and the 
onset of gas accretion onto the cores, we are particularly concerned with the residual Ed- 
distribution at this stage. The asymptotic X^-a distribution is determined by Ci, f g> o and 
Tdep- In this subsection, we show the simulation results with / 9j0 — fd,o, Id = 1-5 and q g = 1.0 
at t — 0. The results with other reasonable values of fd,o/ f g ,o, Qd and q g are qualitatively 
similar. With these initial conditions, we can compute the emergence of cores during the 
epoch of gas depletion. 

Figure [2h shows at t = 10 5 , 10 6 , and 10 7 years (dashed, dotted, and solid lines) with 
Ci = 1 and fg t0 = 3. These results correspond to surviving protoplanets at t ~ r dep for 
Tdep = 10 5 , 10 6 , and 10 7 years, although depletion of f g on time scales Td ep is not taken into 
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account here (in the Monte Carlo simulations in section 3.2, the exponential decay, eq. (J7|), 
is considered). We generated 1000 semi major axes with a log-uniform distribution in the 
range of 0.05-50 AU and simulated the growth and orbital migration of cores there. In 
these simulations, gas accretion onto cores is neglected in order to clearly see the effects 
of type-I migration. The dynamical interac tions among the cores is also neglected. With 



N-body simulations, iKominami et all (120051 ) showed that the dynamical interactions with 



other cores and planetesimals do not change the cores' type-I migrati on speed. The result s 



presented in this panel confirm the finding of N-body simulation by iDaisaka et all (120061 ). 
They also show that is depleted in an inside-out manner. The clearing of planetesimals, 
along the paths of previous generations of cores, limits the growth of the migrating cores 
and the gravitational interactions between them. In order to take into account the effect of 
dynamical clearing in the Monte Carlo simulations (to be presented in next subsection), we 
terminate the growth of cores after their semimajor axes have decreased below Odep,mig given 
by eq. ([27D. 

Figure [2H shows that the total mass of the remaining population of cores is ~ 0.1M®. 
For these retained cores, r mig i > 3r C;acc (r c>acc = M c /M c ), where a factor 3 reflects an actual 

1/3 _____ _____ 

timescale to reach M c because T cacc oc M c . From eqs. (fT3"j) and (TH|) . the maximum mass 
of remaining cores (M c max ) is given by 

m_ * o, lc - (iff (*fcf (^r (0 v <*> 

The above expression with f gfi = f dj0 reproduces the result at > ade P ,mi g ~ 1AU in Figure [2H 
as well as the weak dependence of M c max on f g and C\ (Figures [3bl andHU). In order to 
further examine the / 9i0 an d C\ dependences, we carry out a set of models with / ff)0 = 30 
{Ci = 1) and C\ = 0.1 (/ 9i o = 3) in Figures [3] and HI In the former case, the disks is 
marginally self-gravitating and contain a significant fraction of the central stars' mass. The 
emergence of relatively massive cores is more sensitively determined by C\, the inefficiency of 
type-I migration. A large amount of solid materials in the disk does not efficiently promote 
the formation of massive cores, because the associated dense gas increases type-I migration 



speed, as already pointed out by IDaisaka et all (120061 ) 



In the inner regions at a < 0.3AU, the cores' accretion proceeds on very short time 
scales and they reach their isolation mass near their birth place (Figure [2b). Thereafter, 
most of these cores migrate into their host stars within 10 5 years. Consequently, the local 
Erf is essentially depleted by the formation and migration of the very first generation cores 
(panel c). This domain is determined by the condition M c max > M c iso . From eqs. (TTTT) and 
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f l23|) . this condition implies 

a < ^ = 0.45C- (f) ^ (^) 14/15 AU. (24) 

This expression approximately reproduces the result in Figure [2] and it accounts for the 
dependence on C\ and / 9j o (see Figures [3b andHb). The location of ai SO ,mig does not vary 
with the slope of the surface density distribution because both competing processes are 
determined by local properties of the disk. 

In the limit of efficient type-I migration (with C\ = 1), most of the cores undergo orbital 
decay before they attain their isolation mass in the intermediate regions (for Figures[2j 0.3AU 
< a < 10AU). Consequently a significant fraction of remains to promote the formation 
of subsequent-generation cores. The results of our simulations (Figure [2b) show that many 
generations of cores may emerge at the same disk location. This repeated formation and 
self-destruction process is even more efficient in massive disks where / Sj o ^> 1 (see Figure 
Eb). 

In the inner region, surface density is significantly depleted by the N generations of core 
formation and disruption where 

W " M c , max - ^ \ 3 j \ f g , J \1AU) \M e 

As the planetesimal building blocks become depleted, the cores formed at later epochs have 
masses smaller than M Cjmax . Therefore the above expression for N geaej i slightly under esti- 
mates the number of populations of cores which may emerge. 

At large disk radii, the number of generation is limited by disk depletion time. In this 
region, 

M t 1/4 /4o\ 13/10 fv^Uo \ 3/4 /_a_\ -W ( t \ /AO - 1/4 

W r c , mig (M c , max ) - 6 ^ \ 3j \ f 9 .o J \1AUJ Wye a rs)\M Q 

(26) 

This estimate completely reproduces the results in Figures Eb, Eb, and Hb. The actual 
number of generation is given by min(iVg enei i, iV geil e,2)- Significant depletion of the original 
inventory of heavy elements occurs in the limit N genet x < iV genei 2, that is, within the location 

. < ^ = (if) 1/40 (^) " 4 g) V8 AU. ,27) 

This boundary of disruption zone is in excellent agreement with the critical location within 
which Erf (equivalently, from f^) has reduced from its initial values by an order of magnitude. 
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Note that the dependences of ad ep , m ig on C\ and / Sj o are very weak. As long as iV genei 2 > 1, 
the disruption zone is confined to a ~ (t/10 6 years) 1/ ' 4 AU, independent of disk surface density 
and migration speed C\ (in the limit of small Ci, N genC)2 < 1 and iV genejl < 1, so that no 
depletion occurs). 



These results indicate that, during the early epoch of disk evolution, cores form and 
migrate repeatedly to clear out the residual planetesimals. This self-regulated process pro- 
vides an avenue for the host stars to acquire most of the heavy elements retained by the 
planetesimals. In large, well-mixed molecular clouds where star clusters form, this self- 
regulated mechanism would lead to stellar metallicity homogeneity (§4.1). The inner disk 
region contains cores which have started their type-I migration but not yet reached to the 
disks' inner edge. The average value of at these locations, including the contribution of 
these migrating cores, is two orders of magnitude smaller than its initial value. 

Eventually, the disk gas is so severely depleted that relatively massive cores no longer 
undergo significant amount of type-I migration. As any given disk radius, the condition 
for retaining 90% of the initial solid surface density is N gene> 2 — 0.1iVg enel . We find, from 
eqs. (1251) . (126"j) . and (1271) . that this condition is satisfied in regions with 



O ^ ^survjinig 



ion 



^dep,mig — 2(2 t jep,mig- 

Type-I migration leads to a transition in the distribution at this orbital radius. 



(2? 



Inside the ice line, type-I migration limits the mass of individual surviving cores. But 
these cor es can coalesce through giant impacts durin g and after the severe depletion of the 



disk gas (jKominami fc Ida! 120021 ; llda fc Lin I l2004al ). Provided the total mass of residual 



planets may form near 1AU. Previous simulations f 


Chambers & Wetherill 


1998; 


Agnor et al. 


1999; 


Kominami & Ida 


2002: 


Ravmond et al. 


2004 


) show that the most massive terrestrial 



planets tend to form in inner regions of the computational domain where the isolated cores 
are initially placed. In these simulations, strong gravitational scattering process can inject 
planetesimals close to the host stars to form smaller planets with relatively close-in orbits. 

In general, type-I migration leads to clearing of planetesimals close to their host stars 
and sets the inner edge of the cores' population at ~ 1AU. The lack of pla nets inside the 



Merc ury's orbit in our Solar system might also be attributed to this result (IDaisaka et al. 



20061 ). In addition, the self-regulated clearing process also leads to fd ~ 0(1) near 1AU for 
wide variety of initial conditions (fd,o)- This residual distribution of heavy elements at Td ep 
ensures the formation of Mars to Earth-sized terrestrial planets in habitable zones (see §3.2). 
Even in the limit that the disks' initial distribution is much larger than that of MMSN, 
the reduction of at < a few AU inhibits in situ formation of gas giants interior to the ice 
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line. In contrast to the results in Papers I-III, the inclusion of type-I migration suppresses 
the rapid and prolific formation of hot Jupiter which in term facilitates the formation and 
retention probability of habitable terrestrial planets. 

We also simulated several models with qj = 2 which correspond to a steeply declin- 
ing initial surface density distribution. Such a steep gradient of could be produced 
by the inward migration of dust due to the h ydrodynamic drag on them by the disk gas 



(IStepinski fc Valageaslll997l ; iKornet et al.l 120011 ). In this model, a significant amount of solid 



mass is contained in inner disk regions where cores quickly form and undergo type-I migra- 
tion. Nearly all the initial mass of solid components in disks is accreted by the host stars 
through the self-regulation by type-I migration. The results of this model is consistent with 
the discovery of metallicity homogeneity among the stars in young open clusters (see §4.1). 

The above discussions indicate that the formation of Mars to Earth-sized habitable 
planets depends only weakly on type-I migration speed. However, it is critical for the for- 
mation of cores of gas giant planets, because M Cjmax oc C x 3 ^ 4 (eq. [23]). For giant planets to 
actually form, sufficiently massive cores must be able to accrete gas on time scales at least 
shorter than a few folding time of r dep . For r dep ~ 10 6 -10 7 years, we deduce, from eq. (IToT) 
with kl — 9 and k2 = 3, that gas giant formation is possible only for M c > a few M®. 
But, type-I migration suppresses the emergence of such massive cores in disk regions with 
relatively large E 5 . The results with C\ — 1 in Figures [2H and(3H show that the asymptotic 
core masses are generally much smaller than that needed to launch efficient gas accretion 
even though there is little decline in the magnitude of S^. 

In the case of C\ = 0.1 (Figure HH), the cores formed at a few AU originally have 
~ M©. Provided the planetesimals along their migration path are not captured onto their 
mean motion resonances, these cores can grow up to a few M ffi through accretion during 
migration. The core mass is marginal for rapid gas accretion. As shown in eq. f !23p . in 
metal-rich disk regions (where fd,o/f g ,o > 1), more massive cores can be retained with the 
same amount of solid materials. In this expression, the effects of gas depletion has not 
been taken into account. But, this process is included in the numerical simulations in §3.2 
and it also enhances the mass of retainable cores, especially those which emerge during the 
advanced st ages of gas depletion . These results suggest that gas giants can be formed for 



C\ < 0.1. (jAlibert et al.l (120051 ) derived a similar condition for formation of Jupiter and 
Saturn.) In the next subsection, we will show that C\ < 0.1 reproduces mass and semi 
major axis distribution of extrasolar gas giants comparable to those observed. 
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3.2. Mass - semimajor axis distributions 

3.2.1. Dependence on type- 1 migration speed 

We now consider the signature of type-I migration on the M p — a distribution for a 
population of emerging planets. In the Monte Carlo simulations, we first generate a 1,000 
set of disks with various fd,o (the initial value of fa) and Td op . We adopt the same prescriptions 
for the distributions of fd,o and f g> o as those in Paper II. For the gaseous component, we 
assume f 9j o has a log normal distribution which is centered on the value of / S) o = 1 with 
a dispersion of 1 (<51og 10 / 9)0 = 1.0) and upper cut-off at / Sj0 = 30, independent of the 
stellar metallicity. For the heavy elements, we choose fd,o = 10' Fe / H ' d / 5i o, where [Fe/H]^ is 
metallicity of the disk. We assume these disks have the same metallicity as their host stars. 

Following our previous papers, we also assume r dep has uniform log distributions in 
the ranges of 10 6 -10 7 yrs. For each disk, 15 a's of the protoplanetary seeds are selected 
from a uniform log distribution in the ranges of 0.05-50AU, assuming that averaged orbital 
separation between planets is 0.2 in log scale (the averaged ratio of semimajor axes of adjacent 
planets is ~ 1.6). This procedure is the same as that adopted in Paper II. Constant spacing 
in the log corresponds to the spacing between the cores is proportional to a, which is the 
simplest choice and a natural outcome of dynamical isolation at the end of the oligarchic 
growth. The log constant spacing for planets and cores with similar masses also maximizes 
dynamical stability. In the present paper, we neglect dynamical interaction between planets 
(this issue will be addressed in future papers) and the growth of individual planets are 
integrated independently. Although the choice of averaged orbital separation is arbitrarily, 
it would not change overall results with regard to the effects of type-I migration. 

In all the simulations presented here, a = 10~ 3 and M* = 1M are assumed. Since the 
on-going radial velocity surveys are focusing on relatively metal-rich stars, we present the 
results with [Fe/H] = 0.1 in most cases. The dependence on [Fe/H] is shown in Figures In 
order to directly compare with observations, we determine the fraction of stars with currently 
detectable planets as r]j. In the determination of r]j, we assume the detection limit is set by 
the magnitude of radial velocity (v r > lOm/s) and orbital periods (TV < 4 years). According 
to the following uncertainty, we exclude close-in planets with a < 0.05AU in the evaluation 

of T}j. 

We artificially terminate type-I and II migration near disk inner edge at a 2 day period 
(~ 0.03AU for M* = 1M ) in a similar way to Papers I-III. We have not specified a survival 
criterion for the close-in planets because we do not have adequate knowledge about planets' 
migration and their interaction with their host stars near inner edge of their nascent disks. 
Hence, we record all the planets which have migrated to the vicinity of their host stars. 
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In reality, a lar ge fraction of the gia nt planets migrated to s mall disk radii may either be 
consu med (e.g., Sandquist et al.|[l998l ) or tidally disrupted (e.g., Trilling et al. 1998 ; Gu et al. 



20031 ) by their host sta rs. Cores that have migrated to inner edge may also coagulate and 



form super-earths (e.g.. lTerquem &: Papaloizoull2007l ). We also neglect such core coagulation 
near inner edge. 

In a set of fiducial models, we adopt M* = 1M Q , [Fe/H] = 0.1 and (kl, k2) = (9, 3) for 
the gas giants' growth rate (eq. [T7]). The analytic deductions in the previous subsection 
suggest that relatively massive cores can be retained to form gas giants provided C\ < 0.1. 
Figures [5] show the predicted M p -a distributions for C\ — (panel b), C\ = 0.01 (panel 
c), C\ = 0.03 (panel d), C\ = 0.1 (panel e), and C\ = 0.3 (panel f). In order to directly 
compare the theoretical predictions with the observed data, we plot (in panel a) M p which 
is a factor of 1.27 times the values of M p sini determined from radial velocity measurements 



( |http: / / exoplanet.eu/ ) . This correction factor corresponds to the average value l/(sinz) = 
4/7r for a sample of planetary systems with randomly oriented orbital planes. To compare 
with the theoretical results with M* = 1M , we plot only the data of planets around stars 
with M* = O.8-1.2M observed by radial velocity surveys. 

All results show "planet desert," which is a lack of intermediate-mass (M p ~ 10-IOOMq) 
planets at < a few AU. However, formation probability of gas giants dramatically changes 
with C\. The fraction of stars with gas giants rjj changes from 22.9% for the model with 
C\ = to 0.2% for G\ = 0.3 (see Table [Tj). In the observed data, r\j ~ 5-8% around 
[Fe/H]~ 0.1 (Figure ED- 

In contrast, the distributions of retained gas giants are similar to each other except 
for relatively large C\. The theoretical predictions are also consistent with the observed 
distribution in panel a. We carry out a KS test for statistical similarity between the simulated 
models and the observed data for the parameter domain of 0.1AU < a < 2.5 AU and 
M p > 100M e . This range corresponds to a maximum rectangular region in which planets 
are detectable by radial surveys with precision v r > lOm/s and duration < 4 years. Since 
we have not imposed any criterion for determining the survival probability of short-period 
planets, the predicted population of planets with a < 0.1AU are excluded in the quantitative 
statistical significance test. Except for the model with C\ = 0.3 (panel f), these models 
are statistically similar to the observed data within a significant level Qks ^ 0.3 for both 
semimajor axis and mass cumulative distribution functions. In particular, the model with 
C\ = 0.03 (panel d) shows an excellent agreement with Qks — 0.86 for the mass function. 

In models with C\ = 0.3, only the low-mass cores can survive type-I migration. The 
envelope contraction time scales for these low-mass cores are generally much longer (eq. [IT] ) 
than the gas depletion time scales. Consequently, rjj is very small (< 1%) for the simulated 
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model with C\ = 0.3. Since type-II migration occurs after planets have acquired a mass 
which is adequate to open up gaps, the close-in planets with > 100M e are rare for models 
with C\ = 0.3, in contrast to the models that type-I migration is neglected or sufficiently 
reduced (panels b to e). 

Assuming the survival fraction of close-in planets is independent their M p , the magnitude 
of G\ can be calibrated from the observed mass distribution close-in planets. In Figure [TJi, 
we plot the mass function for all the planets which are halted artificially at a = 0.03 AU. In 
this panel, we neglect any further evolution including both disruption and collisions. We also 
consider an alternative limit that after the gas depletion, the multiple-generation of cores 
which migrated to the proximity of any given star are able to merge into a single terrestrial 
planet with a mass M mer (Fig. [7b). These results suggest the potential findings of many hot 
earths, including Neptune-mass planets, with either transit or radial velocity surveys. 

The above results clearly highlights the competing effect of type-I migration and gas 
accretion. In these models, we approximate the gas accretion process with jkl, k2) = 
(9,3). The early models of proto-gas-giant planet formation (jPollack et al.l Il996l ) yield 
slower growth rates and they are better fitted with (kl,k2) = (10,3.5). But, recent re- 
visions ( Ikoma et al. 2000 ; Hubickyj et ali 2005 ) of these models indicate that the protogas 
giants' growth rates can be significantly enhanced by the opacity r eduction associated with 
grain growth or boundary conditions at different regions of the disk ( llkoma et al.ll200ll ). The 
Kelvin-Helmholtz contraction timescale may also be reduced by turbulent heat transport in 
the outer envelope of the protoplanets. In view of these uncertainties, we also simulated 
models with kl = 8 and 10, with k2 = 3 for all cases. The predicted r]j are listed in Table 
[TJ For models with C% = 0.03-0.1 in which type-I migration marginally suppresses the for- 
mation of gas giants, the magnitude of rjj depends sensitively on the minimum mass for the 
onset of dynamical gas accretion (which is represented by kl). The results in Tabled] indicate 
that a smaller value of kl(= 8) can lead to a significant increase in rjj because smaller mass 
cores can initiate the runaway gas accretion within Td cp = 10 6 -10 7 years (eq. [TT]). 



3.2.2. Preservation of terrestrial planets 



Although type-I migration is an effective mechanism to suppress the emergence of gas 
giant planets, earth-mass planets can form at ~ 1AU even in the limit of G\ ~ 1. (But most 
terrestrial planets with a < 0.3AU would be eliminated by a full-strength type-I migration.) 
Not all of these planets formed interior to the ice line may be retained. A fraction of these 
planets may be trapped onto the m ean motion reson ance of migrating gas giant planets 
and be forced to migrate with them (IZhou et al.ll2005l ). In the relatively unlikely event of 



inefficient eccentricity damping and very fast type-II migration, some of the embryos along 
the path of the migrating gas giants may also be scattered to large radial distances and form 
later generation terrestrial planets ([Raymond et al.ll2006l ). 



The survival of terrestrial planets depends on their post-formation encounter probability 
with migrating giant planets. In the absence of any type-I migration, this probability is 
modest. But the inclusion of a small amount of type-I migration significantly reduce the 
fraction of stars with massive close-in gas giants because the retention of the progenitor 
cores becomes possible only at the late stages of disk evolution when the magnitude of E 9 is 
reduced. With a limited supply of the residual disk gas, the growth of gas giants and their 
type-II migration are suppressed. 

We find that repeated migration of gas giants is less common in models with C\ > 0.01 
than those with C\ = 0. The low type-II migration probability reduced the need for efficient 
disruption of largely accumulated close-in planets (see Paper II). It also ensures that most 
of the terrestrial planets formed in the habitable zones are not removed by the migrating 
gas giants. Note that type-I migration also inhibits in situ formation of gas giants near 1AU 
(see §3.1). Thus, a small amount of type-I migration facilitates formation and retention of 
terrestrial planets in habitable zones in extrasolar planetary systems, rather than inhibits 
them. 



3.2.3. Metallicity dependence 

We also study the dependence on other parameters, metallicity ([Fe/H]). As fiducial 
models, we set C% = 0.03 and kl = 9 to compare the metallicity dependence. The most 
massive cores which can form and be retained prior to gas depletion have masses M Cjmax 
given by eq. fl23l . Smaller / Sj0 , larger /^o, or smaller C\ increases M C)Daax and would enhance 
formation of gas giants. The models in Figures [5] have already illustrated the dependence on 

In Paper II, we showed that relatively large values of [Fe/H] (or equivalent fd,o/ f g ,o) 
enhance the growth rates and increase the isolation masses for the cores, for the disks with 
the same f 9t0 . Fast emergence of massive cores can lead to the rapid onset of gas accretion. 
Based on that model, we found rjj increases with [Fe/H], and the simulated of rjj depen dence 



on the [Fe/H] is qualitatively consistent with observed data (IFischer fc Valenti II2005I ). 



Here we study the effect of type-I migration on the r/j-[Fe/H] correlation. In Figures 
[Fe/H] = 0.1 is assumed. We carried out similar simulations with various [Fe/H] and the 
simulated r/j-fFe/H] relation is plotted in Figure E]). A comparison between the results here 
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and those in Paper II shows that type- 1 migration enhances the ?yj-[Fe/H] correlation. In 
metal-poor disks, fd,o/ f g ,o < 1 an d f]j is significantly reduced by type-I migration. But in 
more metal-rich disks, relatively massive cores can be retained before the disk gas is severely 
depleted. The resultant steep dependence is in a better agreement with the observed data 
(open circles in the plot) than in the dependence without type-I migration (filled circles with 
dashed line), although different choice of assumed averaged orbital separation may change 
r]j slightly. 

Formed planetary systems are affected by stellar mass M* as well as metallicity. The 
dependence on M* was studied in Paper III, using a simple model without the effects of type- 
I migration. In the paper, we predicted that gas giants are much more rare around M type 
dwarfs than around FGK dwarfs while super-Earths are abundant a round M type dwarfs , 



which is consistent with radial velocity survey and microlensing survey (IBeaulieu et al.ll2006l ). 
Adding the type-I migration to the simple prescription in Paper III, we found that the above 
conclusions do not change. Around M type dwarfs, super-earths at 1-3 AU, which are inferred 
to be abundant by microlensing survey, survive type-I migration. In a separate paper, we 
will address the details of the M* dependence of planetary systems, taking into account the 
M* dependences of many physical quantities. 



4. Summary and Discussions 

In our previous Monte Carlo simulations of planet formation processes (Papers I, II, and 
III), we neglected the effects of type-I migration. In the present paper, we have investigated 
its effects on formation of terrestrial planets and cores of gas giants. We found that type- 
I migration provides a self-clearing mechanism for planetesimals in the terrestrial-planet 
region. Although the planetesimal disk at a < 1AU is significantly cleared, the total mass 
of residual planetary embryos at regions within a few AU is comparable to the Earth's 
mass, almost independent of the disk and migration parameters. Earth-like planets can be 
assembled in habitable zones after the depletion of the disk gas. 

But this self-regulated clearing process does prevent giant planets from forming at ~ 
1AU even in very massive disks. In general, the clearing of cores leads to the late formation 
of gas giants. When the surface density of the disk gas is reduced below that of the MMNM, 
type-I migration would no longer be able to remove cores which are sufficiently massive 
to initiate the onset of rapid gas accretion. This late-formation tendency also reduces the 
fraction of gas giant which undergo extensive type-II migration. Since migrating gas giants 
capture and clear cores along their migration paths, type-I migration also facilitates the 
retention of Earth-mass planets in the habitable zones. 
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In the limit that type-I migration operates with an efficiency comparable to that deduced 
from the traditional linear torque analysis (i.e. with C\ — 1), all cores would be cleared prior 
to the gas depletion such that gas giant formation would be effectively suppressed. However, 
the catastrophic peril of type-I migration would be limited by at least a ten-fold reduction 
in the type-I migration speed (C\ < 0.1). In this limit, a substantial fraction of the cores 
may survive and gas giants would form and be retained with an efficiency (rjj) comparable 
to that observed. 

The above discussions indicate that the formation probability of gas giants is delicately 
balanced by various competing processes. Since these processes have comparable efficiency, 
the fraction of solar type stars with gas giants appears to be "threshold" quantity. Small 
variations in the strength of one or more of these effects can strongly modify the detection 
probability of extrasolar planets. For e xample, y.j is observed t o be a rapidly increasing 



function of their host stars' metallicity (IFischer fc Valenti II2005I ). Although we were able 



to reproduce this observed trend in Paper II, the simulations presented here indicate that 
this effect is enhanced by a small amount of type-I migration because it is more effective in 
metal-poor disks. 

The M p - a distribution in Figures [5] also show its sensitive dependence on various other 
model parameters. The most noticeable dependence is the relative frequency between gas 
and ice giant planets at several AU from their host stars. These distributions also predict 
a modest population of close-in cores, a fraction of which may survive and be observable as 
hot earths. 



4.1. Metallicity homogeneity of open clusters 



In modern paradigm of star formation (IShu et al.l 1 198 71 ). most of the stellar content 
is processed through protostellar disks. Gas diffusion in these disks is regulated by the 
process of turbulent transport of angular momentum whereas the flow of heavy elements is 
determined by orbital evolution of their main carriers, i. e. through gas drag of grains and 
tidal interaction of planetesimals and cores with the gas. In general the accretion rates of 
these two components are not expected to match with each other. In fact, the formation of 
gas giant planets around ~10% of the solar type stars and the common existence of debris 
disks requires the retention of heavy elements with masses at least compara ble to that of 



2002; 


Quillen 


2002; 


Shen et al. 


2005) 



The observationally determined upper limit in the 
metallicity dispersion (< 0.03-0.04 dex) among the stars in the Pleiades and IC4665 open 
clusters implies a total residual heavy element mass (including the planets) to be less than 
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twice that in Solar system planets. 

Clues on the resolution of this paradox can be found in the protostellar disks. Recent 
models of the observed millimeter continuum SED's suggest that some fraction of classical 
T Tauri disks has du st mass > 10 times the total heavy-element mass in Solar system 



( iHartmann et al.ll2006l ). Their host stars would acquire significant metallicity dispersion if a 
large fraction of the heavy elements in these disks is retained as terrestrial planets and cores 
while most of their gas components is accreted by the stars. Hence, an efficient process to 
clear heavy elements is required. 

Although, due to gas drag, dust grains can undergo inward migration and be accreted 
onto their host stars, we here assume that the formation of planetesimals is a more efficient 
self-clearing process. In active protostellar disks where £ 9 and are much larger than 
those of MMSN, planetesimals quickly grow into cores which undergo rapid type-I migra- 
tion. Through a series of numerical simulations, we show that most cores formed interior 
to 

^dep,mig ~ 1 AU are lost before the disk gas is severely depleted. On a larger spatial and 
time scales, this process is also effective in transporting most of the original heavy elements 
to their host stars. Provided their progenitor molecule clouds is thoroughly mixed, this 
self-regulated clearing process would ensure stars in young clusters acquire nearly uniform 
metallicity. 



4.2. Diversity of giant planets around solar-type stars 

Among the stars with known gas giants, a large fraction of them show signs of additional 
planets. This special multiplicity function is associated with the formation of gaps around 
the first-born gas giants. Beyond the outer edge of the gap, a positive pressure gradient 
leads the gas to attain a local super Keplerian velocity. Dust particles and planetesimals 
accumulate in this re gion and grow into s ufficiently massive cores to initiate the formation 



of additional planets (IBry den et al. Il2000l ). It is also possible that the planet formation is 
a threshold phenomenon, i.e. the conditions needed to form a multiple-planet systems are 
marginally more stringent than that for the formation of single gas giants. 

The simulations presented here consider the formation probability of individual planets. 
Although the impact of type-I migration on the residual disk and the formation of multiple 
generation cores have been taken into account, we have neglected the impact of gas giant 
formation on the residual disks. Nevertheless, this algorithm can be used to qualitatively 
describe a threshold scenario for the formation of multiple planets. 

The total formation timescale of gas giants is Tf 0Tm ~ T cacc + tkh- In comparison with 
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both type-II migration and the disk depletion timescales, we find that, provided r c>acc < r m i g i, 



Solar-system-like giant planets would form, if Tf orm ~ Td ep , because gas giants would 
not have enough time to undergo extensive type-II migration. 

Eccentric giants would form, if Tf orm < Td ep and r mig 2 > Tf orm . The second condi- 
tion implies that nearby second-generation gas giants are likely to emerge, either 
by spontane ously satisfying the formation conditions or due to an induced forma- 



tion process (jBryden et al. 1 120001 ) . before they undergo any significant type-II migra- 



tion. Orbital instability of closely packed multiple-planet systems would excite their 

eccen tricities ( Rasio fc Fordll996 ; Weidenschilling fc Marzari 1996 ; Lin fc Idalll997 ; Zhou et al 
2007h . 



Close-in giants would form if Tf orm < Td ep , and r mig2 < T{ 
migration before neighboring gas giants form. 



They undergo type-II 



In the second case, if r mig2 < Td cp and orbital instability does not occur on a time scale 
~ Td ep , the multiple planets could be locked into mean-motion resonances during migration. 
According to the above conditions, the regions in which close-in, eccentric, and solar-system- 
like giant planets are likely to form are schematically plotted in the /d,o~ a o plane (ao is the 
initial semimajor axis) in Figure [8], neglecting type-I migration. Bright gray, dark gray, and 
modest gray regions correspond to close-in, eccentric, and solar-system-like giants regions, 
respectively. A more quantitative set of simulations will be presented in a subsequent paper. 
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k 2 


= 10 


21.3 


12.0 


7.9 


2.3 


0. 



Table 1: i]j (in percent) for various C\ and k±. Other parameters are fixed: metallicity 
[Fe/H]= 0.1, stellar mass M* = 1M . 
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Fig. 1. — Viscous evolution of a disk with a = 10 3 . We set S 9 = at r = 0.1 AU and 
10 5 AU. (a) The distributions of E g at t = 0, 10 3 , 10 4 , 10 5 , 10 6 , and 10 7 years (from top to 
bottom), (b) Time evolution of £ s at planet-forming regions r = 0.25, 0.67, 1.8, 4.7, and 12.5 
AU (from top to bottom). 
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Fig. 2. — The evolution of S d due to dynamical sculpting by type I migration in the case of 
Ci = 1 and / fli0 = 3. Here we assume constant f g . (a) The £ d -distributions at t — 10 5 , 10 6 , 
and 10 7 years are expressed by dashed, dotted, and solid lines. The initial distribution is also 
shown by dotted lines, (b) The evolution of /^-distributions. (The meanings of the lines are 
the same as panel a.) (c) The number of generation of protoplanetary cores formed at each 
a. (d) Planet distributions at t — 10 7 years. The dashed line expresses the core isolation 
mass and the scattering limit. 
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Fig. 3. — The same as Figure [2] except for / s>0 = 30 to approximate the early phases of 
massive disk evolution. 
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Fig. 4 - 



The same as Figure [2] except for C\ = 0.1 to approximate weak type I migration. 
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& | AU] a [AU] 

Fig. 5. — Planetary mass and semi major axis distribution. Units of the mass (M p ) and 
semimajor axis (a) are earth mass (M® = Me) and AU. (a) Observational data of extrasolar 
planets (based on data in |http : / / exoplanet . eu/ [ ) , around stars with M* = 0.8-1.2M Q detected 
by the radial velocity surveys. The determined M p sini is multiplied by l/(sini) = 4/7r ~ 
1.27, assuming random orientation of planetary orbital planes, (b) The distribution obtained 
from Monte Carlo simulations without taking into account the effect of type-I migration, 
(c) that includes the type-I migration with C\ = 0.01, (d) C\ = 0.03, (e) G\ = 0.1, and (f) 
C\ = 0.3. The dashed lines express observational limit with radial- velocity measure precision 
of v r = lOm/s. In these models, M* = 1M , the magnitude of the metallicity [Fe/H] = 0.1 
and the contraction time scale parameters in eq. [T7] are assumed to be (kl, k2) = (9, 3). 
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Fig. 6. — The metallicity dependence of fraction of stars with planets which are detectable 
by radial velocity surveys with measurement precision v r > lOm/s and coverage p eriods < 4 
years. Open circles with error bars are observed data (jFischer fc Valenti II2005I ). Dashed, 
solid, and dotted lines with filled circles are the theoretically predicted dependences with 
C\ = 0, 0.03 and 0.3, respectively. Other parameters are (kl, k2) = (9, 3), and M* = 1M Q . 
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Fig. 7. — The predicted mass function which are halted artificially at a = 0.03 AU. (a) 
The mass function for all the planets, neglecting further evolution due to both disruption 
and collisions is neglected. Filled circles, open circles, squares, and triangles represent the 
results with d = 0,0.001,0.01, and 0.1, respectively. (M* = 1M , [Fe/H] = 0.1 and 
(kl,k2) = (9,3).) (b) That includes the maximum effect of coagulation after the gas is 
depleted. All the planets are coagulated into one planet in each system. 
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Fig. 8. — The parameter domains (initial disk surface density fd,o and initial semimajor 
axis of planets a ), in which close-in, eccentric, and solar-system-like giant planets are likely 
to form, in the limit of negligible effects of type-I migration. Bright gray, dark gray, and 
modest gray regions correspond to close-in, eccentric, and solar-system-like giants regions, 
respectively. 



